% Estimate effects of US monetary policy shocks using model from
% Antolin-Diaz and Rubio-Ramirez (2018) under different combinations of
% (narrative) sign restrictions. Goal is to: 1) assess sensitivity of
% results to choice of conditional prior; and 2) disentangle
% informativeness of different restrictions.

clear variables
close all
clc

oldFolder = pwd;
cd ..
addpath([oldFolder,'/auxFunctions']);
cd(oldFolder);

%% Import data.
% Data are obtained from the replication files to Antolin-Diaz and
% Rubio-Ramirez (2018).
% The variables are:
% 1) Real GDP 2) GDP deflator 3) Commodity price index 4) Total reserves
% 5) Nonborrowed reserves 6) Federal funds rate
load('Uhlig_Data_Updated.mat');

% Re-order variables so that the shock of interest is ordered first. 
data = 100*data(:,[end 1:end-1]); % Federal funds rate ordered first

% Declare any exogenous variables (other than the constant; leave empty 
% if none).
exog = [];

% Set variable names for plotting figures.
varnames = {'Federal Funds Rate','Real GDP','GDP Deflator',...
    'Commodity Price Index','Total Reserves','Nonborrowed Reserves'};
% Set units for impulse responses. 0 = %, 1 = ppt, 2 = bps.
irfUnits = [2 0 0 0 0 0];

%% Options.
opt.p = 12; % No. of lags in VAR
opt.const = 1; % const = 1 if constant in VAR, = 0 otherwise
opt.jshock = 1; % Index of shock(s) of interest
opt.ivar = [1 2 3 4 5 6];  % Indices of variables of interest
opt.cumIR = []; % Indices of variables for which cumulative impulse response should be used
opt.H = 60; % Maximum horizon for impulse responses
opt.draws = 1000; % No. of draws from posterior of reduced-form parameters (phi) with non-empty identified set
opt.aalpha = 0.68; % Credibility level for credible intervals (0<aalpha<1)
opt.qMax = 100000; % Max no. of draws of orthonormal matrix Q when checking for empty identified set
opt.qDraws = 1000; % No. of draws of Q used to approximate bounds of (conditional) identified set
opt.dispIter = 100; % Print number of draws of phi remaining every dispIter draws
opt.gridLength = 1000; % Size of grid used when computing credible intervals

rng(23032021); % Set seed for random number generator

restr.elastRestr = []; % No elasticity restrictions

%% Sign restrictions only
% Each row of signRestr contains a vector (i,j,h,s,t) representing a
% `traditional' sign restriction, where t is the type of restriction:
% t = 1: the impulse response of the ith variable to the jth shock at the 
% hth horizon is nonnegative (s = 1) or nonpositive (s = -1).
% t = 2: the (ij)th element of A0 is nonnegative (s = 1) or nonpositive 
% (s = -1). 
% signRestr = []; % No sign restrictions on impulse responses
restr.signRestr = ...
      [1 1 0 1 1; % Response of FFR to monetary policy shock on impact is nonnegative
       1 1 1 1 1; % As above after one month
       1 1 2 1 1; % As above after two months
       1 1 3 1 1; % As above after three months
       1 1 4 1 1; % As above after four months
       1 1 5 1 1; % As above after five months
       3 1 0 -1 1; % Response of GDPDEF to monetary policy shock is nonpositive
       3 1 1 -1 1; % As above after one month, etc
       3 1 2 -1 1;
       3 1 3 -1 1;
       3 1 4 -1 1;
       3 1 5 -1 1;
       4 1 0 -1 1; % Response of COMM to monetary policy shock is nonpositive
       4 1 1 -1 1; % As above after one month, etc
       4 1 2 -1 1;
       4 1 3 -1 1;
       4 1 4 -1 1;
       4 1 5 -1 1;       
       6 1 0 -1 1; % Response of NBR to monetary policy shock is nonpositive
       6 1 1 -1 1; % As above after one month, etc
       6 1 2 -1 1;
       6 1 3 -1 1;
       6 1 4 -1 1;
       6 1 5 -1 1]; 

%% Run again with full set of sign and narrative restrictions.
% Each row of shockSignRestr contains a vector (i,s,t) representing the
% shock-sign restriction that the ith shock in period t is >= 0 (s = 1) or
% <= 0 (s = -1). t is equal to the index of the observation.
restr.shockSignRestr = ...
    [1, 1, find(dates == datenum(1974,4,01));
     1 1 find(dates == datenum(1979,10,01));
     1, 1, find(dates == datenum(1988,12,01));
     1, 1, find(dates == datenum(1994,2,01));
     1, -1, find(dates == datenum(1990,12,01));
     1, -1, find(dates == datenum(1998,10,01));
     1, -1, find(dates == datenum(2001,4,01));
     1, -1, find(dates == datenum(2002,11,01))];

% Each row of hdSignRestr contains a vector (i,j,t,h,k,s) representing a
% restriction on the historical decomposition. 
% k represents the class of restriction (Type A or Type B) and s represents
% the type of restriction within the class.
% k = 1: Type A restriction
%  - s = 1:  The jth shock is the 'most important contributor' to the 
%   change in the ith variable between periods t and  t+h.
%  - s = -1: The jth shock is the 'least important contributor' to the
%   change in the ith variable between periods t and t+h.
% k = 2: Type B restriction
%  - s = 1: The jth shock is the 'overwhelming contributor' to the change
%   in the ith variable between periods t and t+h.
% - s = -1: The jth shock is a 'negligible contributor' to the change in
%   the ith variable between periods t and t+h.
restr.hdSignRestr = ...
    [1, 1, find(dates == datenum(1974,4,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1979,10,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1988,12,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1994,2,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1990,12,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1998,10,01)), 0, 1, 1;
     1, 1, find(dates == datenum(2001,4,01)), 0, 1, 1;
     1, 1, find(dates == datenum(2002,11,01)), 0, 1, 1];
% First shock in these periods was most important contributor to change in
% first variable (federal funds rate).

% Conduct (robust) Bayesian inference.
fprintf('Estimating model with sign restrictions, NR 6 and NR 7');

mainfile;

